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Projection stabilisation applied to general Lagrange multiplier finite element methods is in- 
troduced and analysed in an abstract framework. We then consider some applications of the 
stabilised methods: (i) the weak imposition of boundary conditions, (ii) multiphysics coupling on 
unfitted meshes, (iii) a new interpretation of the classical residual stabilised Lagrange multiplier 
^ , method introduced in H. J. C. Barbosa and T. J. R. Hughes, The finite element method with 

yp ■ Lagrange multipliers on the boundary: circumventing the Babuska-Brezzi condition. Comput. 
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^ ■ I. INTRODUCTION 

O 

<N 

The use of Lagrange multipliers to impose constraints in the finite element method is 
a well-known and powerful technique. To obtain a stable method the finite element 
spaces for the primal variable and the multiplier must be carefully matched so as to 
J^ , satisfy an inf-sup condition uniformly in the mesh parameter (see Babuska [T], Brezzi 

[9], Pitkaranta [30], [31]). If for some reason an unstable pair is used, stability can be 
recovered using a stabilised method [31 135] . 

In many cases such as when imposing incompressibility for flow problems there are 
several choices available, both to design inf-sup stable velocity-pressure pairs (see for 
instance [TU]) and to design stabilised methods for pairs that do not satisfy the inf- 
sup condition. A class of method that has been particularly successful are projection 
stabilisation methods. Loosely speaking such methods ensure stability by adding a term 
that penalises the difference between the pressure solution and its projection onto some 
inf-sup stable space [TT| 16] ITS] [T2] . 

Recently there has been renewed interest in Lagrange multiplier method in the con- 
text of imposing constraints on embedded boundaries and multi-scale or multi-physics 
coupling problems [5] [ITl [U [29l [28] . Also here care must be taken to chose pairs of finite 
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element spaces that satisfy the appropriate inf-sup condition, in order to avoid spurious 
oscihations or locking. 

In some of these cases, although the choice of stable space is known, it may be incon- 
venient. Either the spaces may be very complicated to design or use from an implemen- 
tation point of view, or the multiplier space simply is too small to give sufficient control 
of the constraint. Here the state of the art method for stabilisation is the residual based 
formulation introduced by Barbosa and Hughes jSj. This method has been shown to 
be closely related to Nitsche's method, in cases where the Lagrange multiplier can be 
eliminated locally [33]. It can also be applied for interface coupling proplems, with a 
large flexibility in the choice of multiplier space, see for instance [26 . 

It appears that the idea of projection stabilisation, that has been very successful 
for Stokes' problem, has not yet been exploited to its full potential in the context of 
other type of problems featuring Lagrange multipliers. However it appears that such an 
approach can give certain advantages. 

• For domain decomposition with non-matching meshes it allows for the use of a 
Lagrange multiplier that is defined on a third mesh which can be chosen arbitrarily 
(typically structured). In this case the stabilisation operator only acts on the 
multiplier space, see [H]. This reduces the problem of interpolating between two 
fully unstructured meshes to that of interpolating from two unstructured meshes 
to one structured mesh. 

• Another example is fictitious domain methods where the multiplier can be chosen 
piecewise constant per element and distributed in the interface zone if projection 
stabilisation is used [IS]. This choice is advantageous from the point of view of 
implementation, but normally prohibited since the inf-sup condition fails [23] . 

• Compared to Nitsche type methods or the Barbosa-Hughes stabilised method the 
projection stabilised multiplier method does not use the normal trace of the stress 
tensor explicitly. This is particularly advantageous in the nonlinear case, since the 
nonlinearity then appears only in the bulk and not in the interface terms. 

Stabilised Lagrange methods seem to be attracting increasing attention, in particular 
for the imposition of embedded Dirichlet boundary conditions [TT] [2TJ [TBI ISZl H]. It 
is interesting to note that the extension to XFEM type interface coupling methods is 
practically always straightforward. 

The focus of the present paper is on the generality of this type of method. We 
prove a wellposedness result for discrete solutions and a best approximation result in an 
abstract framework. Then we show how to apply the ideas to the analysis and design of 
stabilised Lagrange multiplicator methods first in the simple case of the weak imposition 
of boundary conditions and then sketching an unfitted finite element method for multi- 
physics coupling. 

As a last example of the applicability of our framework we give a new interpretation 
of the methods of Barbosa & Hughes [3] and Baiges et al. [2], for the imposition of 
boundary conditions. In these methods, the stabilisation acts on the difference between 
the multiplier and the gradient of the primal variable. Using a recent stability result 
for the penalty- free, unsymmetric Nitsche's method [14], we show that such stabilised 
Lagrange multiplier methods are in fact closely related to projection stabilisation methods 
by the inf-sup stability of the Lagrange multiplier space consisting of normal gradients 
of the primal variable on the boundary trace mesh. 
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As a model problem the reader may consider the Poisson problem set on an open 
connected domain ft C M'', d = 2,3, with polygonal (or polyhedral) boundary. Find 
M : J7 — > R such that 

—Au ~ f in n 

(1-1) 
u = on oil. 

The weak formulation of this problem, using Lagrange mulitpliers to impose the boundary 
constraints take the following form: find {u,X) G _ff^(0) x H^^i^Q^) 

Vu • Vf dx + Xv ds + fiu ds = fv dx 

n Jan Jan Jn 

'i(v,n)eH\n)xH-^dn). (1.2) 

We will frequently use the notation a "^b for a < Cb where C is a constant independent 
of the mesh-size, but not necessarily of the local mesh geometry. We also assume quasi- 
uniformity and shape regularity for all meshes. 



II. ABSTRACT SETTING 

We will here give an abstract framework for this type of method to give some under- 
standing of the underlying idea. Our aim is to make the simplest possible framework. 
Let 

a(-,-) :y X T/^M 

and 

6(-,-) -.LxV^M. 

be two bilinear forms representing the partial differential operator on weak form and the 
constraint respectively. The abstract formulation then writes: Find {u,X) £V x L such 
that 

a{u,v) + b(\,v) + b(pL,u) = {f,v) (2.1) 

for all (w,/i) ^ V X L. We assume that the spaces V and L are chosen such that the 
problem is well posed. Firstly we assume that the bilinear forms satisfy the following 
continuities 

a{u,v) < \\u\\v\\v\\v, 

b{\v)<\\M\L\\v\\v 
and secondly that the form a{u, v) is coercive on the kernel of 6(A, v), i.e. 

||f lly < a{v,v), for all v such that b{ji,v) — Q, V/i G L. 

Finally we assume that the Babuska-Brezzi condition is satisfied so that VA G L there 
holds 

I|A||l < sup-^-— . 
vev \\v\\v 
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Example 1. In the case of the Poisson problem (|l.ip above the bilinear forms are given 
by the weak formulation (jl.2p as 



and 



a{u,v) :== / Vu-Vv dx, (2.2) 

Ja 



b{X,v) := / Ai; ds. (2.3) 

Jdn 



The spaces are given by V := H^{n) and L := H 2 (957). 



Now consider the discretisation of the problem (|2.1|) in Vh C V,Lh C L and we assume 
that these spaces satisfy the discrete version of the inf-sup condition uniformly so that 
VXfi G Lfi there holds 

Uhh^ < sup -| — 1^ , (2.4) 

VheVh \\Vh\\Vh 

where || • \\l^ and || • Hv^ are norms defined on the discrete spaces. It is known [251 [^^ that 
the discrete inf-sup condition is equivalent to the existence of an interpolant np : V ^ Vh 
such that \/v ^ V there holds 

b{v -TTpv,^ih) ^0 y^ih ^ Lh, and\\7Tpv\\v <\\v\\v 

Example 2. For V = H^{il) and L = H^^l^Qft) we have 

WtJ-hUn ■= Wh'^l^hWmdn) 
and 

hhWvh ■= l|Vw/i||L2(t^) + \\h^^Vh\\L^dn)- 
It follows immediately by the Cauchy-Schwarz inequality that the following continuities 
hold 

0'iuh,Vh) <\\uh\\vh\\vh\\vh yuh,VheVh (2.5) 

and 

H^h,vh)<\\Xh\\Ljvh\\v, yXheLh^yvhGVh. (2.6) 

This leads to the following formulation: find {uh, Xh} G Vh x Lh such that 

a{uh,Vh) + b{Xh,Vh) -b{uh,Hh) ^ {f,Vh), y{vh,lJ.h} <^Vh y- Lh- (2.7) 

Then we know that the discrete problem is well posed and we may prove optimal con- 
vergence provided the spaces have optimal approximation properties. We will denote the 
kernel of 6(-, •) by 

Kh := {vh e Vh ■■ b{nh, Vh) = 0, V/x/, £ Lh}. 

In the following we will assume that {u, X} are sufficiently regular that ||m||v^ and ||A||i,fe 
are well defined. We will also assume that || • ||y^ is equivalent to || • ||y on V. 

Consider now the case where we do not want to use the space Lh because it leads to 
inconvenient interpolation problems or insufficient stability. We want to work with the 
possibly completely unrelated, richer, space Ah, for which no stability is known to hold, 
but which is convenient from the point of view of implementation. Instead of (|2.4[) we 
then have teh following stability property. 
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Lemma 2.1. For all X^ G A;j there holds 

ll\ II <- H^h.Vh) II I 

IIAhlli^ < sup + \\Xh-nLXh\\L^, 

vhev,, \\Vh\\vH 

where tt^ : A^ -^ Lh denotes the projection operator from A^ to Lh obtained by minimiz- 
ing \\\h -ttlKWlu- 

Proof. By the triangle inequality there holds 

IIAhlli^ < ||A,, - TTiAftllL^ + \\lTLXh\\L^- 

Since 7rz,A?j G L^ the condition (|2.4p holds and hence 

ii\ II ^ ii\ \ II I b{^L>^h,Vh) 
IIA/iIIl;, < IIA/i-ttlA/Jl, + sup —r — 

Vh£Vh \\'"h\\VH 

II \ \ II , Ht^lK ~ Xh,Vh) + b{Xh,Vh) ,^„s 
= ||A/i -vrLA/illL^ + sup — . (2.8) 

VheVh \\Vh\\Vh 

Applying now the continuity of 6(-, •) we conclude 

ii\ II <r ii\ \ H I b{^h,Vh) 
Uh\\L^<\\Xh~TTLXh\\L^+ sup . 

VheVh \\Vh\\Vh 

■ 
This means that provided that we can control the distance \\Xh — T^L^hUbh from the 
approximation in the space Ah to the space Lh, which satisfies the LBB-condition, then 
we will have stability using the space A^ . The simplest way of obtaining this is to add a 
symmetric operator —s{Xh, fJ-k), designed so that 

||A^-7riA/,||i,^<s(A^,A,0 (2.9) 

to the formulation (|2.7p . We also assume that s(-, •) satisfies 

< s{fih,f^h) < CsWtihWl^ for aU ^h G A^. (2.10) 

Here the constant Cg can be made small by chosing the stabilisation parameter in s(-, ■) 
small. Since the effect of s(-, ■) is to reduce the effective dimension of the space Ah it can 
be thought of as a coarsening operator. 
This leads to the stabilised formulation: 

a{uh,Vh) + 6(A/i, Vh) + b{^h,Uh) - s{Xh,^J,h) = if, Vh) 

for all {vh,^ih) G Vh X A,,. (2.11) 

The signs in (|2.1ip have been chosen so as to preserve symmetry, note however that 
the problem is indefinite due to the saddle point structure. For the operator s(-, •), the 
following design criteria are advantageous: 

• minimal dependence of the stable subspace Lh 

• the smallest possible stencil 

• optimal weak consistency. 

Often s{-, ■) may be chosen as the jump of the function or of function derivatives over 
element faces in the multiplier space. Currently this choice, when applicable, is the one 
that to the best of our knowledge, satisfies the above criteria best. 
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When we work with the multipher space A;j, it is no longer sufficient to assume that 
a{ufi,Vfi) is coercive on the kernel of b{\h,Vh), for A/j £ K^. Indeed the stabilisation 
term could upset the coercivity. To ensure that the constraint remains strong enough 
compared to the penalty term so as not to compromise coercivity we assume that for all 
Uh G Vh there exists £,h{uh) G A/i such that 

\\S,h{uh)\\L^ < c^\\uh\\v,, 

(2.12) 
a^lluhlly^ < a{uh,Uh) + b{S,h{uh),Uh). 

S.hi'U'h) is typically related to the constraint that one wishes to impose. For the case of 
weak boundary conditions typically S,h{uh) is the projection of the trace of Uh onto the 
Lagrange multiplier space as we shall see later. 

We conclude this section by the following theorem, showing a best approximation 
property for the formulation (|2.1ip . 

Lemma 2.2. For all {u^jAh} G Vn x Lh there holds 
\\uh\\v^+s{\hAh) < a{uh,Uh)+biXh,Uh)-b{\h-^h{uh),Uh) + s{Xh,Xh-Oi{uh))- (2.13) 
Proof. Starting from the right hand side of (|2.13l) we have using (|2.12p 

a{uh,Uh) +b{Xh,Uh) - b{\h - S.hiuh),Uh) + s{\h,\h - S,h{uh)) 

> oiiW^hWv^ + -s{Xh,Xh) - ■^s{^h{uh),^h{uh))- 
Using now (12.101) followed by the first inequality of (I2.12p . leading to 

we may conclude, assuming c^ small enough. ^ 

Theorem 2.3. Assume that the coercivity condition (I2.12|) holds for Vh x Ah and that 
there exists a space Lh such that the condition (|2.4p holds for the pair Vh ^ Lh- 
The system (j2.1ip admits a unique solution {uh,Xh}- Moreover there holds 

ll'"-'"/i||v^ + 1|A- A,,||l^ < ini \\u~yh\\vn+ inf {\\X - i^hh^ + s{i^h,'^h)'^)- 

Vh&Vh fhGAh 

Proof. Assume that Uh and Xh exist. Now by the triangular inequality 

||u - UhWvn < \\u - T^Fu\\v^ + WtTfU - UhWvh, 

where ttf is the Fortin interpolant associated to the spaces Vh x Lh, and 

||A - A^IIl, < ||A - i^hh,, + \Wh ~ A/,||l^, \lvh G A,,. 
Set r\h — Uh — TTpu and ^h — Xh ~ Vh, by Lemma l2.2l we have 

\\rih\\l,+s{(:h,Ch) 

< a(?7h, Vh) + b{Ch,Tlh) - b{Ch ~ UVh), Vh) + s{Ch, Ch - UVh))- (2.14) 

Subtracting ()2.1ip from (|2.ip with v — Vh, fi = fih, gives the Galerkin orthogonality 

a{u-Uh,Vh) + 6(A- Xh,Vh) + b{^h,u-Uh) + s{Xh,Hh) = 0. (2-15) 
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Taking Vh = ?7/i and /i/j = —{Ch — £.h{Vh)) in (|2.15p and adding the left hand side of (|2.15p 
to the right hand side of (I2.14p yields 

WWiWv^ + siChXh) < a(u - nFU,rih) + b{X - Vh,rih) 

~ KCh - ^hiVh), U - TTpu) + s{vh, Ch - £,h{Vh))- (2.16) 

Since &(/ift, u — t:fu) — for all fih G Lh there holds h{'KLCh, u — ttfu) — and hence 

ll'7/illvfe +s(C/i,C/i) ^ a{u-TrFU,r]h) + b{X - Vh,Vh) + b{^h{T]h),u - irpu) 

- b{Ch ~ ttlCh^u- TTpu) + s(vh,Ch ~ ih{rih))- 

By the continuity of a(-, •) and &(•,•) and the stability of the first equation of (|2.12p . we 
have 

hh\\v^+s{C,hXh) < {\\u-TlFu\\vf, + \\\-Vh\\L^)\\rih\\v^ + \\u~TlFu\\vACh-'^LCh\\Lf, 

+ s{uh,Vh)'^{s{ChXh)'^ + s{^h{w^),ih{llh))'^)■ (2.17) 
Using the upper bound \\C,h ~ T^hChWhu ^ s{ChXh) and 

s{ih{m.),^h{'nh))"^<\\Uvh)\\L,<hh\\v, 

we obtain 

hhfv^ +s{Ch,Ch) < {\\U~TTFU\\V^ + \\\-Vh\\Lu +s{vh,Vh)'^) 

x(||r7/.fy,+s(a,a))*- (2.18) 
This gives the following upper bound for ||?7/i||vfe 

hhWvn +s{Ch-,C,h)^ < \\u-iTFu\\v,, + inf \\\-Vh\\L,, + s{vh,Vh)^ . 
By the stability of ttf we have, for vt e Vu 

\\U - TTFU\\V^ < \\U - Vh\\Vh + \\vh ~ T^FuWvh 

= II" - VhWv^ + hF{vh - u)\\v,^ < \\u - VhWv^. (2.19) 
We conclude that 

\\u-ui,\\v,,+s{ChXh)^< mi \\u-Vk\\v^+ inf {\\X - i^hWl^ + s{iyh,'^h)^)- 

For the bound on (h we use the result of Lemma (j2.ip : 

11/- II <: KCh,Vh) , II . . II ^ b{Ch,Vh) , /^ ^ nI 
\\Ch\\Ln< sup - — + IIU-TTLChlUh ^ sup - — 7 hs(Ch,O0'- 

Vh^Vh \\Vh\\VH VHfEVh \\Vh\\VH 

Since we already have the desired bound for the stabilisation term we only need to 
consider the first term. By the Galerkin orthogonality p.isp . with /i/j = we have 

b{Ch,Vh) = b{\ - Vh.Vh) + a{u - Uh,Vh) < (||A- Vh\\Ln + \\u - Uh\\vJ\\vh\\vh- 
This leads to the upper bound on HC/illih: 

iiaiu. < ||A - i^hh, + \\u- uhWv, + s(a, a)i (2.20) 
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Combining this with the previous bounds on ||u — w;j||y,^ and s{(^h,Ch)^ concludes the 
best approximation resuh. 

To prove the existence and uniqueness of the discrete solution Uh and A^ set / = 
and then prove that Uh — and A^ = 0, implying that the system matrix is regular. 
Since the continuous problem (|2.ip is well posed m = and A = 0, but then by choosing 
yh = and Vh — we have 

WuhWvf, + Uhhi, < inf ||w,,||y,^ + inf \\i'h\\L,,=0, 



Remark. Using the continuous inf-sup condition it is easy to show that convergence 
in the discrete norm |j • ||i^ implies convergence in the continuous norm |j • \\l. Following 
\31f . using the continuous inf-sup condition, Galerkin orthogonality and continuity, we 
have 

b{\-\h,v) b{X- Xh.v -ihv) - a{u-Uh,ihv) 

U - ^hh < sup -— = sup 

vev \\v\\v vev \\v\\v 

< ||A- Aftlli,^ + \\u-Uh\\v^. 

Here ih '■ V -^ Vh ClV denotes a V-stable quasi interpolant with optimal approximation 
properties. 

Remark. Note that the convergence order of the method will depend on the approxi- 
mation properties of the finite element spaces and the weak consistency of the stabilisation 
operator. If this operator can be designed so that s(/i/j,A/i) = for /i^ G L^ then the 
influence of the stabilisation operator vanishes. 



III. PROJECTION PENALTY OPERATORS 

The design of the stabilisation operator s(-, ■) is important, indeed if the construction of 
the projection ttl requires a too detailed understanding of the inf-sup stable space Lh 
the advantages of the stabilised method may be lost. Typically this is the case if ttl is 
chosen to be the L^-projection, which is the natural choice. Fortunately there are some 
operators that can handle a relatively large set of spaces {L;i,A/i}. The two natural 
choices are local projection stabilisation or interior penalty stabilisation. Herein we will 
only discuss the second choice. For examples of local projection stabilisation methods 
that can be used in this context we refer to |14| [7] where such methods have been proposed 
in a different context. The extension to the present case is straightforward. Below we 
will focus on the construction relevant for weak imposition of boundary conditions. The 
extension to domain decomposition is straightforward. 

We consider only one side F of the polygonal boundary dVt. Denote the trace mesh of 
Vh by Fy. Let the space Lh be defined on a trace mesh F^, 

Lh :- {lh e C"(Fi) -.IhW^ Pk{K), VK e F^}, 

and A/,, on a trace mesh Fa, 

A,, := {A,, e L^Ta) : Xhlx e Pk{K), VK e Fa}. 
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The mesh function on the trace meshes wih be denoted hpx, with X = V, L or A. We 
assume that there are positive constants Ci and C2 such that 

^r,y(2;) < Ci/ir,L(x) < C2h-p^A{x), for all a; e F. 

For imposition of boundary conditions and more generally for domain decomposition 
methods the classical condition for inf-sup stability is that /ir,L > Chry for some con- 
stant C > 1, (see jlj). Using the projection stabilisation this condition may be relaxed for 
the space Ah, since the stabilisation controls the unstable modes. The relative difference 
in mesh size should be accounted for in the stabilisation parameter to tune the constant 
of the inf-sup condition. Numerical evidence however indicate that this dependence is 
relatively weak. Assume for simplicity that /ir.A < /ir,L- Let I^ ■ A/i — >■ Lh denote a 
quasi interpolation operator such that for Uh G Ah 

lLUh{xi) :^ N^'^ ^ Uh{xi)\K, for aU nodes Xi of F^, 
{KerA-.xeK} 

where A^^: denotes the cardinality of the set {K : x £ K}. For A/j G A/j it follows that 

Now consider any element in the trace mesh Ky G F^ and map it to a reference element 
Kyj^ ■ Also map the same subset of Fa so that Fa := {K' : K' n A'r ^0}- It then follows 
that 

k 

\\\H-iL\h\\L,^i,,< J2 T.[ KXhTds, 

where [x] denotes the jump of the quantity x over an element face, with [x] = for faces 
on the boundary and 9^ denotes the normal derivative of order n, with 9JJ, defined to be 
the identity. 

This upper bound on the reference element follows by the observation that if the jump 
of Xh and all its normal derivatives are zero, then Xh is a polynomial over all of Kr, but 
since ttl interpolates this polynomial {Xh — ttlA/i)!^ =0. After scaling back to physical 
space and summing over all elements in F^ we obtain 

The order sq depends on Lh and follows from the scaling argument, if the Lh norm is 
equivalent to the L^-norm over F then so = 1. We conclude that the interior penalty 
stabilisation operator may be written 



s{Xh,^Jih) :=J2[ h'"^'^[dlXh]K^ih]ds 



It may be inconvenient to compute all normal derivatives up to polynomial order 
and an equivalent local projection approach may be used instead as suggested in the 
references given above. A minor modification of the above argument allows for spaces 
Lh including discontinuous functions. 
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IV. PENALTY STABILISATION OF LAGRANGE MULTIPLIER FORMULATIONS: 
APPLICATIONS 

As an example of the above theory we recah a method introduced as a fictitious domain 
method in [16\ Here we will first present the method in the simple case of weak imposition 
of boundary condition and then propose an extension to unfitted finite element methods. 
Both cases are considered in two space dimensions, but the extension to three dimensions 
is straightforward. 

A. Weak imposition of boundary conditions 

In this section we will consider the problem (|l.ip . that we recall here for the readers 
convenience. 

Let n be a bounded domain in R^, with convex polygonal boundary dfl. The Poisson 
equation that we propose as a model problem is given by 

-^" = ^ ^" ^l, (4.1) 

u ^ g on oil, 

where di^ denotes the boundary of the domain fl, f € L'^{Vt) and g G H^{drt). Under 
these assumptions (|4.ip has a unique solution u G H^{V.) r\H'^{ft) satisfying ||u||_ff2(o) < 

ll/IU=(o)- 

The usual L^-scalar product on the domain 51 will be denoted by {■,-)n or on the 
boundary < •, • >ao. We also introduce the discrete norms 



and 






Recall that there holds 



{\^^)an<\\M\-^M.^n\\^^\\lh.^n■ (4-2) 

The weak formulation of the problem is given by (|1.2p . Let a(-, •) be defined by (j2.2p and 
&(•,•) by dZl- 



Finite element formulation We introduce a triangulation 7^, fitted to the boundary of fl. 
The set of faces of triangles that form the boundary dil of H. is denoted T. 

We will use the following notation for mesh related quantities. Let hx be the diameter 
of K and h = mayiKeTh ^K- We introduce the finite element spaces 

Vh := {v e H^in) -.vIkG Pi{n), 'iK e %} 

and 

A^ := {^l e L\dn) -. h\f g PoIf), vf g t}. 



It is known that this choice of spaces does not satisfy (|2.4p . 

The standard finite element formulation writes: find u^ G Vh and A^ G A;i such that 

a{uh,Vh) + b{Xh,Vh) + b{ph,Uh) = {f,Vh) + b{^ih,g) for all (vh.Hh) G 14 x Ku. (4.3) 
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Assume that L^ denotes a coarsened version of A^, L}^ C A;j such that the inf-sup 
condition is uniformly satisfied for the pair V;j x i/j, we now that this is always possible 
if Lfi is chosen small enough. As proposed in the previous section we may stabilise the 
formulation (|4.3p by adding the penalty term 



S{\h: IJLh) = {h{\h - TTlA/i), ^Xh " T^L^-h) 



on 



Assume that L/j is the set of piecewise constant functions on a boundary mesh with 
mesh-size H > 3h. This space is known to satisfy the inf-sup condition, see [53]. Clearly 
the space A/j is more convenient to work in since it does not require any special meshing 
of the boundary. If we now let X := {xi} be the set of all the points in dfl excluding 
corner nodes. Then there holds, by the arguments of Section UlI. I 

\\Xh-7r^M\\,^^^,,< J2 h'[Xh]\%. 

XjGX 

This prompts the stabilisation operator 

Xj<EX 

Leading to the stabilised method: find Uh € Vh and A;i G A/j such that 

a{uh,Vh) + b{Xh,Vh) + b{nh,Uh) -js{Xh,^ih) 

= (/, Vh) + b{nh, g) for ah {vh, fJ-h) e V^ x A^. (4.4) 

We will outline the analysis of the penalty stabilised Lagrange multiplier method using 
the abstract framework derived in Section IlITI 



Satisfaction of the assumptions of the abstract analysis We may now use the abstract 
analysis of Theorem 12.31 combined with Lemma |2. II to prove a best approximation result 
using the norms 

WuhWv,, = \\uh\\iji: ||A/J|l^ :== ||A||_i^;,_ao- 

By assumption Lh satisfies the inf-sup condition (|2.4p and hence we have the stabilised 
inf-sup condition (12. ip . It is easy to see that the continuities (|2.5I) and (|2.6I) hold. The 
condition ()2.12p also holds by taking £,h{uh) '■— Sh~^^TroUh, where S € M"*" and TTQUh G A^ 
denotes the facewise projection of Uh on the piecewise constants, defined by 

{nouh, yujon = {uh, yh)on : ^Vh e A,,. 

The satisfaction of (|2.12p now follows from the stability of the L^-projection and the 
definition of the Lh norm, 

\\h-^^GUh\\L^ < \\uh\\i,h^9n < \Wh\\v,- (4.5) 

The second relation of (|2.12p is satisfied using the approximation property of the projec- 
tion TTg 

\\uh - TroUhh^h.,d!n < co||Vm/i X nan\\_ih^QQ 
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and a trace inequality ||Vu,i x n9n||-i,/i,ao ^ ct||Vu/i||L2(s2), leading to 

hhfi^hMn ^ W^^h-TTQUhWlj^g^ + WnoUhW^^i^Q^ < ||7roM/i|||_^gf2-cgc?||Vw,,||i2(o)- (4.6) 
It follows that 

and hence for S < Cq a[ the coercivity assumption is satisfied. 

We conclude that the assumptions of Theorem l2.3l are satisfied and that the formula- 
tion (J4.4I) is wellposed and satisfies a best approximation result. 

B. Unfitted finite element methods and multi-model coupling 

Here we will consider the coupling of two models of elasticity over an interface that is 
not fitted to the computational mesh. This type of method can be useful for problems 
where the interface itself is an unknown and repeated computations have to be performed 
with different interface positions, for instance for transient problems where an interface 
moves through the mesh or for inverse identification where the interface will move during 
iterations. 

We consider a geometrical setting with two subdomains, ili and D,2 and a separating 
interface T . In each subdomain Vii we consider the following partial differential equation: 

V • (Tiiui) = /, in fl, 

where ui £ [H^{fli)]'^ denotes a displacement field, ai{ui) G [iJ(div; ili)]^ the stress tensor 
and / G L^(J7) the applied force. Across the interface we assume that the following 
matching conditions hold 

Ml - "2 = 0, (cri(wi) - (T2iu2)) • nr = 0. 

For simplicity we assume that m = on the outer boundary dil. Let u £ V :— [iJo(J7)]^ 
denote the function such that u\{n. = Ui, i = 1,2. We assume that the following coercivity 
and continuity holds for the continuous problem. 

2 

a\\ufy <J2{<^^{'u^),^u^)n,, VwGF, (4.7) 



|^(a,(u,),Vw,)oJ <M\\u\\v\\v\\v, Vu,veV. 

If in addition we assume that the ai{ui) are linear in ui, then existence and uniqueness 
of the problem follows from the Lax-Milgram lemma. 



Finite element formulation This time let Th :— {K} be a triangulation of il that is 
constructed without fitting the element nodes or sides to the interface F. We now extract 
two subtriangulations, Ti '■— {K GTh ■ KDfli j^ 0}, i=l,2. We define two finite element 
spaces, one for J7i and one for £72 by 

V, := {v G [H\%)r ■■ v\k G [Pi{K)]^ and «|aonr, = 0}. 
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Let Gh := {K <E Th '■ K nr ^ 0}. On the set G/j we define the foUowmg multiplier space 

A^ := {A^ e [L^{Gh)]^ ■■ Xh\K e [Po{K)]^ VK G G^. 

The finite element method once again is on the generic form find u^ :— {uhi,Uh2} G 
Vi X V2 —: Vh and A^. e A^ such that 

a{uh, Vh) + b{Xh,Vh) + b{fj.h, Uh) - s{Xh, Hh) = (/, Vh) 

ioT all {vh, tJ.h) e Vh X Ah, (4.8) 

where this time 

2 

a{uh,Vh) := ^(a'.i(Mw), Vwj)n,, b{Xh,Vh) = (A^, (wi - W2)}p . 

Since we know that if we instead looked for Xh in a space Lh defined on macro elements of 
diameter > 3/i, the inf sup condition would be satisfied, we get the abstract stabilisation 
operator 

s(A/i, fXh) := {h{Xh - ttlA/i), fih - T^L^^hJY ■ (4-9) 

In practice, since we do not wat to be concerned with the construction of Lh we apply 
the ideas of section IIII.I and instead work with the operator 

s{Xh,lih):= Y] I h'[Xh]-[iih]As, (4.10) 

where [x] denotes the jump of the quantity x over the interface T. 

Remark. Note that although the operator of (|4.9p is defined on T the operator (|4.10p 
is defined on the interior faces of elements in Gh ■ This convenient trick introduced in JlSf . 
allows us to use the mesh structure for stabilisation and we never need to worry about the 
actual intersections of T with element boundaries. Uniformity of the stabilisation relies 
on the mesh regularity. 



Satisfaction of the assumptions of the abstract analysis Similarly as in the case of weak 
boundary condition we introduce the following norms on the discrete spaces 

WMllhx - (/^"'^' ^)r ' ll^ll-i^r = ('^^' ^)r 

and 

2 

Mlh := J2 ll^"»'»lli^(n) + IKoM|||,;,,r- 

2 = 1 

To prove that the hypothesis of Theorem l2.3l are satisfied we chose the norms || • [|y^ and 
II • ||l,^ as follows, 

ll^lk. := hill,/., I|A||l. := ||A||_i,,,r- 
We assume that the following discrete equivalent of (|4.7|) holds for this norm 

2 
ah\\u\\jj^<^{cri{uhi),Vuhi)n. + \\Muh]\\'l^h,r Vu^ £ T4. (4.11) 

i=l 
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The condition (|2.12l) also holds by taking £,h{uh) '■— Sh^^TTQlui^]. By this choice 

Uh{uh)\\L^ = \\h^^Muh]\\-ih,r ^ IkoKlll i,h,r ^ \\uh\\v^ 
and 

2 

a{uh,Uh) +b{^h{uh),Uh) = ^(CTi(u/ii), Vui,i) + |l7ro[u;i]|l| ,^ p > a/i||M,J|y^. 
We conclude that the results of Theorem 12.31 hold in this case as well. 



Remark. Note that we can only control the projection of the jump onto the piecewise 
constant functions in this case, since the trace inequality used in the previous example 
depends on the cut of the interface and the mesh and S might have to be taken very small 
in the analysis indicating a local loss of control of [uh] ■ The projection on the other hand 
can be controlled independently of the mesh interface intersection. Some care must also 
taken when verifying the assumption (|4.1ip . since this step involves a Korn's inequality 
if T is a closed boundary. 

Remark. By using suitable extensions of the solution following {16f and [25}j optimal 
convergence may be obtained for smooth solutions. The conditioning of the system how- 
ever depends on how the interface cuts the mesh and must be handled either following 
the ideas introduced in JT^ or by preconditioning. 

C. Nitsche's method and stabilised Lagrange multiplier methods: a different approach 

The close relation between the residual based stabilised methods for Lagrange multipliers 
as introduced by Barbosa and Hughes and Nitsche's method was discussed by Stenberg 
in [33]. The idea of that paper was that if the Lagrange multiplier can be eliminated 
locally by static condensation Nitsche's method is recovered. Other authors have recently 
discussed the need of penalty for Nitsche's method and its close relation to Lagrange 
multiplier methods, see for instance |T31 [HI [H] . 

Herein we will show the connection between the projection stabilised methods dis- 
cussed above and the residual based stabilisation of the Lagrange multiplier. Let us first 
recall (one incarnation of) the method of Barbosa & Hughes: find {uh,Xh} G V/i x A/j 
such that 

AbhUuIi, >^h), {'"h, IJ'h)] 

:= a{uh,Vh) +b{Xh,Vh) - b{uh,^h) + j {h{\h +Vuh ■ n),^h + Vw/,. • n)g^ 

^{f,Vh),y{vh,fih}eVhxAh, (4.12) 

with a(-, •) and b{-, •) are defined by (|2.2p and (|2.3p . corresponding to the weak impo- 
sition of boundary conditions. Recalling that formally the Lagrange multiplier is given 
by the diffusive flux A — — Vw • n, we immediately conclude that the method is consis- 
tent. Stability is then typically proven by testing with Vh — Uh and /i^ = Xh using the 
positivity of the form to obtain control of ||/i2 A;j||af2 by absorbing all the other terms in 
the stabilisation using the _ff ^-seminorm of Uh over the domain. Control of Uh on the 
boundary is then obtained in a second step by choosing /i/i suitably. 

We will now consider the stabilisation used in (|4.12p as a penalty on the distance to a 
stable subspace. This would mean using the space Nh of normal derivatives of Vh on the 
trace mesh as multiplier space, together with Vh for the primal variabel. Since in that 
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case Vfi and N^ no longer can be chosen independently this method may be written: find 
Ufi e Vfi such that 

ANit{uh, Vh) := a{uh, Vh) + h{-Vuh ■ n, Vh) - b{uh, -Vvh ■ n) 

^if,Vh),yvheVh. (4.13) 

We have eliminated the Lagrange multiplier in the formulation using its equivalence with 
the diffusive flux. Writing out this variational formulation leads to 



Vuh ■ Vu/i dx — Vuh ■ nvh ds + Vf/i • riuh ds — fv dx, 
o Jon Jan Jn 

which we identify as the non-symmetric version of Nitsche's method, without penalty. For 
the argument to make sense we now need a stability result for this method. The question 
of the inf-sup stability of the non-symmetric version of Nitsche's method, without penalty, 
was recently treated in [IS, where the following stability result was proven: 

Lemma 4.1. Let Vh be the standard space of piecewise polynomial continuous finite 
element functions. Assume that every side of the polygonal D, is meshed with at least 
five elements, then for some C ^ Cq > 0, with cq independent of h, but not of the mesh 
geometry, there holds 

II II ^ ANit{uh,Vh) 

\\uh\\i,h< sup — - — , 

VheVh \\Vh\\i,h 

where 

IWhWlh := W^UhWh^^^ + CWh-^UhWh^i)^), C > 0. 

It follows that we have the required stability and we may prove stability of the residual 
based stabilisation using the techniques discussed above. 

Remark. The above lemma can be rewritten as 3wfi e Vh such that 

Cw\\uh\\lM< ANitiuh,Wh) (4-14) 

with Wh :— Uh + CVd, Cw > and 

W^dWiJi <co\\uh\\^^h,dn- (4-15) 

The function ipg ensures the control of the boundary contribution. 



We now give an alternative proof of the equivalent of Lemma 14.11 for the formulation 
(I4.12P using the framework of penalty on the distance to the stable subspace. The result 
holds for multiplier spaces satisfying the following compatibility assumption. 

Assumption [Al]: the following continuity holds for the spaces Vh and Ah. For every 
Vh G Vh there exists Zh{vh) € Ah such that 

biuh^'^Vh-n+Zhivh)) < Cz\\Vuh\\L^n)Wyh\\mn), \\h^Zh{vh)\\L'^(dn) < Cz\\^Vh\\L^n)- 

(4.16) 
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Theorem 4.2. Let Vh x A^ satisfy assumption [Al]. Then for all {u^, A/j} G V^ x A^ 
there holds 

\\uh\\iM+\\Xh\\L^< sup — ^ . — -| — II (4.17) 

Proof. First take Vf^ = u^ and ii^ = A/i to obtain 

= a[uh,Uh) + b{\h,Uh) - b{Xh, Uh) + 7 {K^h + ^Uh ■ n),Xh + Vuh ■ n)g^ (4.18) 
We add and subtract Vu^ • n and Vw/i • n in tlie 6(-, •) forms of the formulation 

a{uh,Vh) + b{\h,Vh) -b(uh,fih) + 1 {HXh + Vu,, ■n),Hh + Vv,i ■n)g^ 

= a{uh, Vh) + b{-Wuh ■ n, Vh) + b{Vuh ■n + Xh,Vh) 

- b{uh, -Vvh ■ n) - b{uh, Vvt ■ n + fih) + 7 (H^ + Vu^ ■ n), Hh + Vw^ • n)gj^ 

= ANit{uh,Vh) 

+ b{Vuh ■ n + A/i, Vh) - Huh; Vvh ■ n + /i/i) + 7 {h{Xh + Vuh ■ n), fih + ^Vh ■ n)g^^ . 

(4.19) 

We will first show that by taking Vh — Wh (of (|4.14p - (|4.15p ) and /i/i := A/j + Zh(Cfd) we 
have 

WuhWuh + iHXh + ^Uh ■ n||i2(ao) < ABH[{uh, Xh), {wh, Xh + ZhiCvd))]- 
First note that by the construction of Wh and /i^ and the form (|4.19p we have 

ABHiiuh, Xh), {wh, t^h)] = ANit{uh,Wh) + j\\h^{Xh + Vuh)|||2(an) 

+ C&(Vu/i • n + A/i, ipd) - C,b{uh,Vipd ■ n + Zh{fd)) 

+ 7C {h{\/uh ■ n + Xh), Vipd ■ n + z/i(¥'a))aa • (4.20) 

Then note that by the continuity of 5(-, ■) and by using the Cauchy-Schwarz inequality 
in the penalty term we have 

Cb{Vuh ■n + Xh,ipd) < Cb\\h^{Xh + Vuh)\\L-2(di-i)C,Cd\\uh\\^^hMn 

< \7\\hH^h + yuh)\\l2^9a)+^''C'4cl\\uh\\lh,an^ (4.21) 

C,b{uh,^LpQ -n + Zhiifd)) < CCz\Wuh\\L^(n)\Wd\\i,h 

< \cA^Uh\\h^n) + Clc-^cie\\nh\\Xh,on (4-22) 
and 

j((h(Vuh ■n + Xh),VipQ ■n + Zh{(pd)) < t7||^"(-^/i + V"'!)llL2(an) 

+ 2{cyg + cl)e7\\uh\\l^^g^. (4.23) 

If we choose C smaU enough, it follows from (|4.20p . (|4.14p - (|4.15p and the bounds (|4.2ip - 
(|T^ that 



ABH[{uh,Xh)Awh,lJ.h)] > xCt„||M/J|?,ft + ^TII^K^/i + V"/i^"^ 



L2(ao)- 
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Since \\h^Xh\\L^dn) < \\h^Xh + '^Uh)\\L^{dn) + ctW'^UhWL^n} we deduce that 

\\uh\\lji + II^^A,i|||2(aj^) < ABH[{uh, Xh), {wh, fJ.h)]- 
It only remains to show that 

\\uh + C'PdWiM- + \\h^{Xh + Zh{C'fid))\\mdn) ^ \\uh\\iji + Wh^XhWmon)- 
This is immediate by the triangle inequality and the stability 

\\^d\\i,h + \\hhhiC(pd))\\mdn) < \\uh\\^^h,dn- 



Remark. The condition (I4.16P is easily satisfied for any reasonable space A/j. For 
spaces including discontinuous functions on boundary elements Fj take Zh{vh)\F- '■= 
Jp Vuh-n ds. If the spaces Af,, consists of continuous functions decompose the boundary 
in macro patches Fj consisting of a sufficient number of elements for the construction of 
functions Zh{vh) G Hq(Fj) such that Jp Zh{vh) ds — Jp Vu/,. ■ n ds. 



Remark. It follows from Theorem \4-S\ that the stabilisation of Barbosa and Hughes, 
can be interpreted as a penalty on the distance to the stable subspace, where the stable 
subspace is the normal derivatives of the primal finite element space in the setting of the 
non-symmetric Nitsche method. 

Remark. Note that formally we only require 7 > 0, however if ^ is chosen too small 
both the control of the Lagrange multiplier and the control of the boundary integral of Uh 
are lost. When 7 is chosen too big on the other hand the control of the boundary integral 
of Uh is lost. This can be avoided if the Lagrange multiplier space is chosen so that 
Zh{'~Pd) -^ ^V(/9a ■ n, i.e. if the space of normal gradients on the boundary is a subspace 
of the multiplier space. Indeed if this is not the case the method can lock for very large 
values of J. 
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